function idat = fsb_normalize_zscore(idat)

% FSB - DEV: Z-transform image data 
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% EXAMPLE:
% idat = fsb_normalize_zscore(idat,sandbox)
%
% INPUT:
% idat:         4-D image data
% sandbox:      sandbox experiment struct
%
% OUTPUT:
% idat:         z-transformed 4D image data
%
% CALLED BY:
% FSB.m
%
% NOTES:
% data transformation is done by subtracting the mean voxel value from the
% individual voxel value and then dividing the difference by the standard
% deviation
%
%$ Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Calculating z-scores...');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Transform int16 to single for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

dim = size(idat);
idat = single(idat);
idat(idat==0)=NaN;
waitbar(10/100,h);
imean = int16(nanmean(idat,4));
waitbar(20/100,h);
istd = int16(nanstd(idat,0,4));
waitbar(30/100,h);

try
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % This is the code path for machines with enough memory
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    imean4d = repmat(single(imean),[1 1 1 dim(4)]);
    waitbar(40/100,h);
    istd4d = repmat(single(istd),[1 1 1 dim(4)]);
    waitbar(50/100,h);
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Do actual z-score transform
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    idat = (single(idat)-imean4d)./istd4d;
    waitbar(90/100,h);
    disp('Using fast code path')
catch
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % This is the slower code path for smaller machines
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    clear imean4d istd4d;
    pack;
    disp('Using alternate code path, this might take a little longer')
    for x = 1:dim(4)
        waitbar(x/dim(4),h);
        idat(:,:,:,x) = (single(idat(:,:,:,x))-single(imean))./single(istd);
    end
    
end

close(h)

end
